rm(list = ls()) # Borramos variables de environment
# descomentar si es la primera vez (y requieren instalación)
# install.packages("tidyverse")
# install.packages("performance")
# install.packages("olsrr")
# install.packages("corrr")
# install.packages("corrplot")
# install.packages("skimr")
library(tidyverse)
library(performance)
library(olsrr)
library(corrr)
library(corrplot)
library(skimr)Entrega I RESUELTA
Instrucciones (leer antes de empezar)
Modifica de la cabecera del documento
.qmdtus datos personales (nombre y DNI). IMPORTANTE: no modifiques nada más de la cabecera.Asegúrate, ANTES de seguir editando el documento, que el archivo
.qmdse renderiza correctamente y se genera el.htmlcorrespondiente en tu carpeta local de tu ordenador (pulsa el botónRendero guarda el documento conRender on Saveactivado).Los chunks (cajas de código) creados están o vacíos o incompletos, de ahí que la mayoría tengan la opción
#| eval: false. Una vez que edites lo que consideres, debes ir cambiando cada chunck a#| eval: true(o quitarlo directamente) para que se ejecuten.Recuerda que puedes ejecutar chunk a chunk con el botón play o ejecutar todos los chunk hasta uno dado (con el botón a la izquierda del anterior).
IMPORTANTE: solo se podrá subir un archivo
.htmlal campus, no se evaluarán entregas sin el.htmlcorrectamente generado. Recuerda: el código es una herramienta, no el fin en esta asignatura. Se evaluará especialmente las interpretaciones y conclusiones que detalles. Un ejercicio con el código perfecto pero sin ningún tipo de razonamiento, interpretación o conclusión no superará el 4 sobre 10 de nota.
Paquetes necesarios
Necesitaremos los siguientes paquetes (haz play en el chunk para que se carguen):
Caso práctico I: análisis de datos de cáncer
El archivo de datos a usar lo cargaremos desde el csv cancer.csv.
# no cambies código
datos <- read_csv(file = "./cancer.csv")
datos# A tibble: 3,047 × 13
deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge Geography
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 165. 61898 260131 11.2 500. 39.3 Kitsap C…
2 161. 48127 43269 18.6 23.1 33 Kittitas…
3 175. 49348 21026 14.6 47.6 45 Klickita…
4 195. 44243 75882 17.1 343. 42.8 Lewis Co…
5 144. 49955 10321 12.5 0 48.3 Lincoln …
6 176 52313 61023 15.6 180. 45.4 Mason Co…
7 176. 37782 41516 23.2 0 42.6 Okanogan…
8 184. 40189 20848 17.8 0 51.7 Pacific …
9 190. 42579 13088 22.3 0 49.3 Pend Ore…
10 178. 60397 843954 13.1 428. 35.8 Pierce C…
# ℹ 3,037 more rows
# ℹ 6 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
# PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
# BirthRate <dbl>
Los datos representan diferentes características socioeconómicas de distintas regiones, extraídas de the American Community Survey (census.gov), clinicaltrials.gov, y cancer.gov.
¿Nuestro objetivo? Ser capaces de predecir la variable deathRate, nuestra variable objetivo, que representa la mortalidad media de cancer, por cada 100 000 habitantes. La idea es que seamos capaz de determinar cual de las variables es la que más efecto (lineal) tiene en la mortalidad de cáncer, realizando el ajuste posterior.
El resto de variables son:
medianIncome: mediana de los ingresos de la región.popEst2015: población de la regiónpovertyPercent: porcentaje de población en situación de pobreza.studyPerCap: ensayos clínicos relacionados por el cáncer realizados por cada 100 000 habitantes.MedianAge: edad mediana.Geography: nombre de la región.AvgHouseholdSize: tamaño medio de los hogares.PercentMarried: porcentaje de habitantes casados.PctHS25_Over: porcentaje de residentes por encima de los 25 años con (máximo) título de bachillerato.PctUnemployed16_Over: porcentaje de residentes de más de 16 años en paro.PctPrivateCoverage: porcentaje de residentes con seguro de salud privado.BirthRate: tasa de natalidad.
IMPORTANTE: siempre que puedas, relaciona los valores numéricos de la salida de R con su fórmula.
Ejercicio 1 (1/10)
Ejecuta el código que consideres para chequear que a) todas las variables son numéricas; b) no hay datos ausentes; c) no hay problemas de rango/codificación. En caso de que alguna de ellas no se cumpla, ejecuta el código que consideres para corregirlo y/o eliminar las variables que consideres. Sé conciso en el código.
Simplemente había que darse cuenta de que
- No había ausentes (no se hace nada)
Geographyera no numérica –> fuera (de momento).- los rangos permitidos:
deathRate: entre 0 e infinito.medianIncome: entre 0 e infinito.popEst2015: entre 0 y 8.000.000.000 (habitantes del mundo xd)povertyPercent: entre 0 y 100 (al ser porcentaje)studyPerCap: entre 0 e infinitoMedianAge: entre 0 y 130 años (por poner un tope)AvgHouseholdSize: entre 0 e infinito (si existe alguna finca enorme)PercentMarried: entre 0 y 100 (al ser porcentaje)PctHS25_Over: entre 0 y 100 (al ser porcentaje)PctUnemployed16_Over: entre 0 y 100 (al ser porcentaje)ctPrivateCoverage: entre 0 y 100 (al ser porcentaje)BirthRate: entre 0 e infinito.
Lo comprobamos con skim() o summary()
datos |> skim()| Name | datos |
| Number of rows | 3047 |
| Number of columns | 13 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 12 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| Geography | 0 | 1 | 16 | 42 | 0 | 3047 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| deathRate | 0 | 1 | 178.66 | 27.75 | 59.70 | 161.20 | 178.10 | 195.20 | 362.80 | ▁▇▆▁▁ |
| medIncome | 0 | 1 | 47063.28 | 12040.09 | 22640.00 | 38882.50 | 45207.00 | 52492.00 | 125635.00 | ▇▇▁▁▁ |
| popEst2015 | 0 | 1 | 102637.37 | 329059.22 | 827.00 | 11684.00 | 26643.00 | 68671.00 | 10170292.00 | ▇▁▁▁▁ |
| povertyPercent | 0 | 1 | 16.88 | 6.41 | 3.20 | 12.15 | 15.90 | 20.40 | 47.40 | ▃▇▃▁▁ |
| studyPerCap | 0 | 1 | 155.40 | 529.63 | 0.00 | 0.00 | 0.00 | 83.65 | 9762.31 | ▇▁▁▁▁ |
| MedianAge | 0 | 1 | 45.27 | 45.30 | 22.30 | 37.70 | 41.00 | 44.00 | 624.00 | ▇▁▁▁▁ |
| AvgHouseholdSize | 0 | 1 | 2.48 | 0.43 | 0.02 | 2.37 | 2.50 | 2.63 | 3.97 | ▁▁▃▇▁ |
| PercentMarried | 0 | 1 | 51.77 | 6.90 | 23.10 | 47.75 | 52.40 | 56.40 | 72.50 | ▁▂▇▇▁ |
| PctHS25_Over | 0 | 1 | 34.80 | 7.03 | 7.50 | 30.40 | 35.30 | 39.65 | 54.80 | ▁▂▇▇▁ |
| PctUnemployed16_Over | 0 | 1 | 7.85 | 3.45 | 0.40 | 5.50 | 7.60 | 9.70 | 29.40 | ▅▇▁▁▁ |
| PctPrivateCoverage | 0 | 1 | 64.35 | 10.65 | 22.30 | 57.20 | 65.10 | 72.10 | 92.30 | ▁▂▇▇▂ |
| BirthRate | 0 | 1 | 5.64 | 1.99 | 0.00 | 4.52 | 5.38 | 6.49 | 21.33 | ▂▇▁▁▁ |
La única con problemas es MedianAge que tiene de máximo 600 años: filtramos por ese umbral de 130 años y los eliminamos (más adelante los pasaremos a NA y los depuraremos).
datos_preproc <-
datos |>
select(where(is.numeric)) |>
filter(MedianAge < 130)
datos_preproc# A tibble: 3,017 × 12
deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 165. 61898 260131 11.2 500. 39.3
2 161. 48127 43269 18.6 23.1 33
3 175. 49348 21026 14.6 47.6 45
4 195. 44243 75882 17.1 343. 42.8
5 144. 49955 10321 12.5 0 48.3
6 176 52313 61023 15.6 180. 45.4
7 176. 37782 41516 23.2 0 42.6
8 184. 40189 20848 17.8 0 51.7
9 190. 42579 13088 22.3 0 49.3
10 178. 60397 843954 13.1 428. 35.8
# ℹ 3,007 more rows
# ℹ 6 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
# PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
# BirthRate <dbl>
Ejercicio 2 (0.8/10)
Realiza un primer análisis exploratorio para entender nuestros datos (cada variable por separado). Hazlo tanto de manera numérica como de manera visual (no hace falta que ejecutas 19 000 gráficos, decide uno de los aprendidos para ello y saca conclusiones, al grano). ¿Hay datos atípicos?
Por ejemplo densidades…
datos_tidy <-
datos_preproc |>
pivot_longer(cols = everything(), names_to = "variable", values_to = "values")
ggplot(datos_tidy) +
geom_density(aes(x = values, color = variable, fill = variable)) +
facet_wrap(~variable, scale = "free") +
theme_minimal()…o boxplots
ggplot(datos_tidy) +
geom_boxplot(aes(x = values, color = variable, fill = variable)) +
facet_wrap(~variable, scale = "free") +
theme_minimal()Si se observan atípicos, en especial en popEst2015 (hay algunas regiones muy grandes), en studyPerCap (es decir, no en todas las regiones se investiga por igual, hay algunas donde se hacen muchos estudios), y las variables de renta (medIncome y PovertyPercent), que muestran una gran desigualdad.
De momento suficiente con darse cuenta de esto, ya intervendremos.
Ejercicio 3 (1.2/10)
Nuestro objetivo será predecir la mortalidad de cáncer. De momento trabajaremos en el contexto de la regresión univariante, así que lo primero que haremos tras preprocesar y entender los datos será seleccionar la variable predictora que mayor efecto lineal tenga respecto a la variable objetivo. Hazlo tanto de manera numérica como visual, de manera que puedas descartar correlaciones espúreas (que no haya un dinosaurio, por ejemplo).
Sacamos primero la matriz de correlaciones
datos_preproc |> correlate()Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'
# A tibble: 12 × 13
term deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 deathRate NA -0.428 -0.119 0.429 -0.0225 -0.00429
2 medIncome -0.428 NA 0.235 -0.788 0.0444 -0.117
3 popEst20… -0.119 0.235 NA -0.0649 0.0557 -0.176
4 povertyP… 0.429 -0.788 -0.0649 NA -0.0562 -0.194
5 studyPer… -0.0225 0.0444 0.0557 -0.0562 NA -0.0317
6 MedianAge -0.00429 -0.117 -0.176 -0.194 -0.0317 NA
7 AvgHouse… -0.0381 0.113 0.109 0.0736 -0.00388 -0.357
8 PercentM… -0.267 0.354 -0.161 -0.642 -0.0383 0.430
9 PctHS25_… 0.405 -0.471 -0.312 0.194 -0.0852 0.331
10 PctUnemp… 0.380 -0.452 0.0519 0.654 -0.0314 -0.128
11 PctPriva… -0.385 0.724 0.0525 -0.822 0.0934 0.0692
12 BirthRate -0.0876 -0.0109 -0.0573 -0.0123 0.0107 -0.101
# ℹ 6 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
# PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
# BirthRate <dbl>
Con focus() del paquete {corrr} podemos filtrar solo la que nos interesa: predictoras vs objetivo
datos_preproc |>
correlate() |>
corrr::focus(deathRate)Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'
# A tibble: 11 × 2
term deathRate
<chr> <dbl>
1 medIncome -0.428
2 popEst2015 -0.119
3 povertyPercent 0.429
4 studyPerCap -0.0225
5 MedianAge -0.00429
6 AvgHouseholdSize -0.0381
7 PercentMarried -0.267
8 PctHS25_Over 0.405
9 PctUnemployed16_Over 0.380
10 PctPrivateCoverage -0.385
11 BirthRate -0.0876
Y podemos visualizarlas
datos_preproc |>
cor() |>
corrplot(method = "color")La variable con una mayor correlación es povertyPercent (en positivo, a más pobreza, más mortalidad) y la variable medIncome (en negativo, a más ingresos, menos mortalidad). El signo importa y su interpretación también.
Antes de quedarnos con una de ellas debemos comprobar que no son correlaciones espúreas, así que visualizamos todas vs objetivo
datos_tidy <-
datos_preproc |>
pivot_longer(cols = -deathRate, names_to = "variable", values_to = "values")
ggplot(datos_tidy, aes(x = values, y = deathRate, color = variable)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap(~variable, scale = "free") +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
Aunque sea una relación débil, no parece aberrante que podamos ajustar una recta a esa nube de puntos (no hay por ejemplo un…dinosaurio). Nos quedaremos por tanto con povertyPercent como predictora
Ejercicio 4 (1/10)
Una vez elegida la variable con mayor efecto lineal, el objetivo será usar dicha covariable para predecir la variable objetivo. Lo que haremos será separar las observaciones en dos datasets distintos: un dataset train con el que entrenaremos el modelo (el que usará la regresión para aprender y determinar sus coeficientes) y un segundo dataset test que usaremos en el futuro para obtener predicciones. Para ello realiza un muestreo aleatorio de manera que train contenga el 80% de los datos y test el 20% restante. IMPORTANTE: no cambies la semilla.
Primero vamos a crear un id de la fila, para luego poder filtrar correctamente. Para ello usaremos rowid_to_column() que nos numera automáticamente las filas.
Después usaremos slice_sample(prop = 0.8, replace = FALSE). Importante: replace = FALSE viene ya así por defecto, pero es importante que entendamos por qué: no quieroe que haya dos observaciones repetidas.
Y por último, para el otro 20% no puedo hacer otro slice_sample ya que no me garantiza que no caigan observaciones que ya están en train. Así que usaremos un anti_join o un filter()
set.seed(12345)
# Creamos id
datos_preproc <-
datos_preproc |>
rowid_to_column(var = "id")
train <- datos_preproc |> slice_sample(prop = 0.8)
# versión antijoin
test <- datos_preproc |> anti_join(train, by = "id")
test# A tibble: 604 × 13
id deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 5 144. 49955 10321 12.5 0 48.3
2 25 237. 37625 23372 22.1 0 41.9
3 39 207. 42121 68714 13.9 160. 42
4 49 210. 35799 27037 20.6 1147. 43
5 53 156. 43835 104236 22.5 95.9 30
6 58 140. 38204 7229 16.7 0 49.1
7 64 169. 37468 29126 20.8 0 43.2
8 85 158. 68430 49762 5.9 0 39
9 103 138. 55955 5937 9.2 0 38.8
10 109 163. 40979 3625 12.7 0 47.5
# ℹ 594 more rows
# ℹ 6 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
# PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
# BirthRate <dbl>
# version filter
test <- datos_preproc |> filter(!(id %in% train$id))
test# A tibble: 604 × 13
id deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 5 144. 49955 10321 12.5 0 48.3
2 25 237. 37625 23372 22.1 0 41.9
3 39 207. 42121 68714 13.9 160. 42
4 49 210. 35799 27037 20.6 1147. 43
5 53 156. 43835 104236 22.5 95.9 30
6 58 140. 38204 7229 16.7 0 49.1
7 64 169. 37468 29126 20.8 0 43.2
8 85 158. 68430 49762 5.9 0 39
9 103 138. 55955 5937 9.2 0 38.8
10 109 163. 40979 3625 12.7 0 47.5
# ℹ 594 more rows
# ℹ 6 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
# PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
# BirthRate <dbl>
Ejercicio 5 (1.8/10)
No usaremos el dataset test hasta el final para evaluar las predicciones (con un dataset que el modelo no conoce). Con el dataset train a) ejecuta el ajuste de una regresión lineal univariante; b) interpreta los coeficientes; c) realiza la diagnosis de la manera más completa posible; d) comenta todo lo que sepas sobre la segunda tabla de la salida (la de inferencia de los parámetros). Elimina variables si su efecto no fuese significativo.
Tras separar los datos, entrenamos el modelo con train
ajuste <- lm(data = train, formula = deathRate ~ povertyPercent)
ajuste |> summary()
Call:
lm(formula = deathRate ~ povertyPercent, data = train)
Residuals:
Min 1Q Median 3Q Max
-121.490 -14.128 1.459 15.520 170.393
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 147.53777 1.44389 102.18 <2e-16 ***
povertyPercent 1.84645 0.07942 23.25 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 25.2 on 2411 degrees of freedom
Multiple R-squared: 0.1831, Adjusted R-squared: 0.1828
F-statistic: 540.5 on 1 and 2411 DF, p-value: < 2.2e-16
De momento lo único que podemos interpretar es
- \(\beta_0 = 147.53777\), es decir, la predicción del modelo da un ratio de mortalidad de \(147.53777\) (por cada 100 000 habitantes) cuando no hay un solo habitante pobre en la región.
- \(\beta_1 = 1.84645\), es decir, por cada punto que suba el porcentaje de pobreza, la mortalidad media por cada 100 000 sube \(1.84645\) puntos.
Así el ajuste formulado sería
\[Mortalidad = 147.53777 + 1.84645 * \%pobreza\]
¿Se cumplen las hipótesis?
check_model(ajuste)A priori parece que gráficamente se cumplen todas. Vamos a chequear una a una
- Errores incorrelados: parece existir independencia entre los errores
check_autocorrelation(ajuste)OK: Residuals appear to be independent and not autocorrelated (p = 0.228).
- Normalidad:
ols_test_normality(ajuste)Warning in ks.test.default(y, "pnorm", mean(y), sd(y)): ties should not be
present for the Kolmogorov-Smirnov test
-----------------------------------------------
Test Statistic pvalue
-----------------------------------------------
Shapiro-Wilk 0.9834 0.0000
Kolmogorov-Smirnov 0.0409 6e-04
Cramer-von Mises 194.4758 0.0000
Anderson-Darling 7.7665 0.0000
-----------------------------------------------
Según el p-valor no hay normalidad…pero vamos a asumir que sí por el siguiente gráfico
ggplot(tibble("residuals" = ajuste$residuals)) +
stat_qq(aes(sample = residuals)) +
stat_qq_line(aes(sample = residuals)) +
theme_minimal()Los valores se pegan a la línea, sobre todo en el centro que es lo importante.
Además si pintamos su densidad vs una densidad de una normal (con misma media y desviación)
set.seed(12345)
ggplot(tibble("residuals" = ajuste$residuals,
"normal" =
rnorm(n = length(ajuste$residuals), mean = mean(ajuste$residuals),
sd = sd(ajuste$residuals)))) +
geom_density(aes(x = residuals), color = "#822181",
linewidth = 1.5) +
geom_density(aes(x = normal), color = "#EFE981",
linewidth = 1.5) +
theme_minimal()Vemos que por lo que puede estar fallando la normalidad es que, por tener poco tamaño muestral E IMPORTANTE NO ESTAR TRATANDO OUTLIERS, la distribución de los errores nos sale ligeramente más apuntada hacía arriba (leptocúrtica, kurtosis positiva).
Esto no había que hacerlo de momento, pero para entender porque los p-valores no concuerdan con lo que ven nuestros ojos, vamos a detectar outliers en la variable predictora
Para eso vamos a usar rápido la función scores() del paquete {outliers}: si ponemos en type = "iqr" (es decir, nos vamos a basar en un boxplot para detectar outliers), nos devuelve cuanto se aleja cada dato de los límites del rango intercuartílico
distancia_x <-
train$povertyPercent |>
outliers::scores(type = "iqr")
distancia_x[1:20] [1] 0.00000000 0.00000000 2.60975610 0.01219512 0.19512195 0.00000000
[7] 0.53658537 0.00000000 1.09756098 0.93902439 -0.37804878 0.00000000
[13] -0.53658537 0.00000000 0.00000000 0.43902439 0.00000000 0.00000000
[19] 0.00000000 0.00000000
distancia_y <-
train$deathRate |>
outliers::scores(type = "iqr")
distancia_y[1:20] [1] 0.00000000 0.08430233 1.51453488 0.71802326 0.55232558 0.00000000
[7] 0.00000000 -0.28779070 0.00000000 0.60755814 -0.08139535 0.00000000
[13] -0.90988372 0.00000000 0.61627907 0.00000000 -1.18313953 -1.70058140
[19] 0.00000000 -0.10465116
Un forma típica y rápida de detectar outliers es cuando dicha distancia supera el 1.5 (en valor absoluto), así que lo que hacemos es eliminar dichas observaciones
train_sin_outliers <-
train |>
filter(abs(distancia_x) < 1.5 & abs(distancia_y) < 1.5)
ajuste_outliers <-
lm(data = train_sin_outliers, formula = deathRate ~ povertyPercent)
ols_test_normality(ajuste_outliers)Warning in ks.test.default(y, "pnorm", mean(y), sd(y)): ties should not be
present for the Kolmogorov-Smirnov test
-----------------------------------------------
Test Statistic pvalue
-----------------------------------------------
Shapiro-Wilk 0.9952 0.0000
Kolmogorov-Smirnov 0.027 0.0676
Cramer-von Mises 186.1388 0.0000
Anderson-Darling 2.941 0.0000
-----------------------------------------------
Empieza a ser normal… Con un tamaño muestral más grande y un tratamiento más rigurosos de los atípicos lo tendríamos.
- Homocedasticidad
check_heteroscedasticity(ajuste)Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).
check_heteroscedasticity(ajuste_outliers)Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).
De la misma manera no parece que se cumpla a priori la homocedasticidad (por lo mismo) pero tenemos lo que nos interesa: errores en torno a una banda constante de varianza
ggplot(tibble("id" = 1:length(ajuste$residuals),
"residuals" = ajuste$residuals),
aes(x = id, y = residuals)) +
geom_point() +
geom_smooth(method = "lm") +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
ggplot(tibble("id" = 1:length(ajuste_outliers$residuals),
"residuals" = ajuste_outliers$residuals),
aes(x = id, y = residuals)) +
geom_point() +
geom_smooth(method = "lm") +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
En ambos casos están en torno a una banda, en el segundo aún más ya que se ha eliminado uno de los outliers en los residuales. Si siguiésemos con un tratamiento mejor de los mismos, voilá.
- Linealidad:
ajuste_res <-
lm(data = tibble("residuals" = ajuste$residuals,
"fitted" = ajuste$fitted.values),
formula = residuals ~ fitted)
ajuste_res |> summary()
Call:
lm(formula = residuals ~ fitted, data = tibble(residuals = ajuste$residuals,
fitted = ajuste$fitted.values))
Residuals:
Min 1Q Median 3Q Max
-121.490 -14.128 1.459 15.520 170.393
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.394e-14 7.713e+00 0 1
fitted -1.395e-16 4.301e-02 0 1
Residual standard error: 25.2 on 2411 degrees of freedom
Multiple R-squared: 2.366e-32, Adjusted R-squared: -0.0004148
F-statistic: 5.705e-29 on 1 and 2411 DF, p-value: 1
ajuste_res <-
lm(data = tibble("residuals" = ajuste$residuals,
"fitted" = ajuste$fitted.values),
formula = residuals ~ fitted + I(fitted^2))
chisq.test(ajuste$residuals, ajuste$fitted.values,
simulate.p.value = TRUE, B = 300)
Pearson's Chi-squared test with simulated p-value (based on 300
replicates)
data: ajuste$residuals and ajuste$fitted.values
X-squared = 784225, df = NA, p-value = 0.003322
Respecto a la linealidad ídem: aunque el p-valor no salga como debería, no se aprecia patrón (ni lineal ni cuadrático)
ggplot(tibble("residuals" = ajuste$residuals,
"fitted" = ajuste$fitted.values),
aes(x = fitted, y = residuals)) +
geom_point() +
geom_smooth(method = "lm") +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
ggplot(tibble("residuals" = ajuste$residuals,
"fitted" = ajuste$fitted.values^2),
aes(x = fitted, y = residuals)) +
geom_point() +
geom_smooth(method = "lm") +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
De hecho fíjate que esa distribución leptocúrtica de nuestros errores hace que nuestra densidad de \(y\) sea un poco más apuntada que las simulaciones (que lo que hacen es usar el ajuste y añadirle el error simulando que fuese normal)
check_predictions(ajuste)Warning: Maximum value of original data is not included in the
replicated data.
Model may not capture the variation of the data.
Fíjate como cambia al tratar mínimamente outliers
check_predictions(ajuste_outliers)Ejercicio 6 (1/10)
¿Se cumplen las hipótesis? Argumenta porqué, más allá de lo que valga luego la bondad de ajuste, es importante que se cumplan.
dsadsa
Respecto a por qué es importante que se cumplan las hipótesis, podemos dar muchos motivos pero principalmente dos vistos en clase con ejemplos simulados:
No cumplir hipótesis implica no poder hacer inferencia: la recta solo nos sirve de manera muestral, para esta no muestra, no pudiendo establecer conclusiones del ajuste a la población. Esto, entre otras cosas, hace que no podamos interpretar ninguno de los p-valores y contrastes, ni que el modelo sirva más allá de esta muestra.
No cumplir hipótesis puede implicar que aunque la bondad de ajuste parezca buena, las predicciones sean bastante malas (por ejemplo, por heterocedasticidad)
Ejercicio 7 (1.4/10)
Evalua el modelo realizando un ANOVA, detallando y explicando todo lo que puedas cada uno de los elementos. Analiza la bondad de ajuste. ¿Qué % de información explica nuestro modelo? Detalla que implica respecto a la calidad del ajuste (en caso de que lo haga).
Para la significación gloval del modelo y el análisis de la varianza usamos anova()
anova_ajuste <- ajuste |> anova()
anova_ajusteAnalysis of Variance Table
Response: deathRate
Df Sum Sq Mean Sq F value Pr(>F)
povertyPercent 1 343293 343293 540.47 < 2.2e-16 ***
Residuals 2411 1531400 635
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
De aquí sacamos
- que la suma de desviaciones al cuadrado de la variable estimada \(\hat{y}\) respecto a su media es de 343293, lo que hemos llamado \(SSR = 343293\)
- que la suma de los residuos al cuadrado es \(SSE = 1531400\)
- que \(SSE / (n-p-1)\) es aprox \(635\), que es lo que hemos llamado como \(\hat{\sigma}_{\varepsilon}^{2}\) (el estimador insesgado de la varianza residual), lo cual cuadra con la anterior salida donde Residual std error era de 25.2, que al cuadrado da 635.04
- además en el contraste global \(H_0:\beta_1=\ldots=\beta_p=0\) (en este caso…\(p=1\)), obtenemos un p-valor muy pequeño –> rechazamos la hipótesis nula, el modelo tiene alguna predictora con un efecto lineal significativo respecot a la objetivo (en este caso, la única que tenemos)
Para el \(R^2\) podemos hacerlo de tres formas:
- con
summary(), que vemos que \(R^2 = 0.1831\)
ajuste |> summary()
Call:
lm(formula = deathRate ~ povertyPercent, data = train)
Residuals:
Min 1Q Median 3Q Max
-121.490 -14.128 1.459 15.520 170.393
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 147.53777 1.44389 102.18 <2e-16 ***
povertyPercent 1.84645 0.07942 23.25 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 25.2 on 2411 degrees of freedom
Multiple R-squared: 0.1831, Adjusted R-squared: 0.1828
F-statistic: 540.5 on 1 and 2411 DF, p-value: < 2.2e-16
- con
compare_performance()aunque solo tengamos un modelo
compare_performance(ajuste, ajuste)Some of the nested models seem to be identical
# Comparison of Model Performance Indices
Name | Model | AIC (weights) | AICc (weights) | BIC (weights) | R2 | R2 (adj.) | RMSE | Sigma
----------------------------------------------------------------------------------------------------------
ajuste | lm | 22425.0 (>.999) | 22425.1 (>.999) | 22442.4 (>.999) | 0.183 | 0.183 | 25.192 | 25.203
- de manera manual con el
anova(), dividiendo \(SSR\) entre \(SST = SSR + SSE\)
anova_ajuste$`Sum Sq`[1] /
(anova_ajuste$`Sum Sq`[1] + anova_ajuste$`Sum Sq`[2])[1] 0.1831195
El \(R^2\) nos dice que nuestro modelo solo consigue explicar el 18% aproximadamente, lo cual puede parecer poco si queremos predecir, pero a nivel médico, implica que solo con tener la variable pobreza de una persona podemos explicar el 18% de su mortalidad de cáncer a futuro (solo con una variable, en un tema tan complejo como el cáncer, not bad)
Ejercicio 8 (0.8/10)
Por último, utiliza el dataset test para, aplicando el modelo entrenado, predecir las observaciones de test. Con las observaciones de test construye un nuevo dataset de tres columnas: la variable objetivo, la variable predictora y la predicción. Haz lo mismo con train. Junta ambos datasets con una cuarta variable que nos diga si la observación era de train o de test.
Lo que hacemos primero es construir un dataset con la variable objetivo y predictora en test, y las predicciones usando predict() (ahora el newdata es literal el conjunto test que hemos apartado).
pred_test <-
tibble("x" = test$povertyPercent,
"y" = test$deathRate,
"y_est" = predict(ajuste, newdata = test),
"split" = "test")Hacemos lo mismo con train y los juntamos (fíjate que he añadido esa cuarta variable para distinguir si las filas vienen de test o de train)
pred_train <-
tibble("x" = train$povertyPercent,
"y" = train$deathRate,
"y_est" = ajuste$fitted.values,
"split" = "train")
pred_data <-
bind_rows(pred_train, pred_test)
pred_data# A tibble: 3,017 × 4
x y y_est split
<dbl> <dbl> <dbl> <chr>
1 13.5 194. 172. train
2 16.1 198. 177. train
3 41.9 248. 225. train
4 20.6 220. 186. train
5 22.1 214. 188. train
6 20 174 184. train
7 24.9 191. 194. train
8 15.7 151. 177. train
9 29.5 186. 202. train
10 28.2 216. 200. train
# ℹ 3,007 more rows
Ejercicio 9 (1/10)
Con el dataset anterior, visualiza de alguna manera la calidad de las predicciones en train y en test (idea: reales vs predicciones) en dos gráficos en el mismo ggplot (pista: la cuarta columna que has añadido antes la puedes usar en un facet). ¿En cuál ha acertado más el modelo? Calcula el \(R^2\) en el dataset de test (piensa como haciendo uso de los errores) y compáralo con el obtenido en la salida del
lm()en train.
Respecto al R2 cuadrado podemos compararlos a mano sabiendo que es siempre ratio de varianza explicada.
R2_data <-
pred_data |>
summarise(R2 = var(y_est)/var(y), .by = "split")
R2_data# A tibble: 2 × 2
split R2
<chr> <dbl>
1 train 0.183
2 test 0.179
IMPORTANTE: no podemos hacer un segundo ajuste en test, y calcular el R2. Lo importante de test es que vamos a evaluar la predicción con el ajuste ya estimado de train, ¡no uno nuevo! queremos “simular” que pasaría si ese modelo ya construido y entrenado, se lo pasamos a unos datos que el modelo nunca conoció
De ahí que salga algo peor en test que en train (lógicamente).
Para el gráfico, vamos a visualizar predicciones vs valores reales, en cada uno
ggplot(pred_data, aes(x = y, y = y_est)) +
geom_point(aes(color = split), size = 2, alpha = 0.7) +
geom_smooth(method = "lm") +
facet_wrap(~split) +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
En ambos casos obtenemos algo que sabíamos: como modelo predictivo, es un modelo de poca capacidad predictiva.
¿Y si metiésemos más variables? Continuará…